Descriptive Information

Load initial script

source("0_load_data.R")

Table describing each variable

# Function to map variable prefix to study wave
get_study_wave = function(var_name) {
  prefix = substr(var_name, 1, 1)
  wave_map = c(
    "a" = "1st Contact",
    "b" = "2 Year",
    "c" = "3 Year",
    "d" = "4 Year",
    "e" = "In Home",
    "g" = "7 Year",
    "h" = "8 Year",
    "i" = "9 Year",
    "j" = "10 Year",
    "l" = "12 Year",
    "n" = "14 Year",
    "p" = "16 Year",
    "r" = "18 Year",
    "u" = "21 Year",
    "z" = "26 Year"
  )
  return(wave_map[prefix])
}

# Function to clean description labels
clean_description = function(x) {
  x %>%
    # Remove range info at end (e.g., ", 0-22", ", 1-9", ", 4-11")
    str_remove(", \\d+-\\d+$") %>%
    # Remove "see value labels"
    str_remove(", see value labels") %>%
    str_remove("see value labels") %>%
    # Spell out abbreviations inline
    str_replace("SOC employment", "Standard Occupational Classification employment") %>%
    str_replace("MHQ\\)", "Mental Health Questionnaire)") %>%
    str_replace("SLQ\\)", "questionnaire)") %>%
    str_replace("G-game", "G-game (general cognitive ability)") %>%
    str_replace_all(" qnr", " questionnaire") %>%
    # Clean up extra whitespace

    str_squish()
}

v_rq1x = 
  data.frame(
    `Short Label` = rq1x_labels_clean,
    Description = rq1x_labels,
    `Teds Code` = ifelse(rq1x %in% colnames(df0), rq1x,paste0(rq1x,"*")),
    `Range or Level` = sapply(rq1x, function(var) {
      if (class(df[[var]]) == "numeric") {
        paste0(round(min(df[[var]], na.rm = TRUE), 2), " — ", round(max(df[[var]], na.rm = TRUE), 2))
      } else if (is.factor(df[[var]])) {
        factor_levels = levels(df[[var]])
        paste(c(paste0(factor_levels[1],"*"), factor_levels[-1]), collapse = ", ")
      } else {
        paste(unique(df[[var]]), collapse = ", ")
      }
    }),
    N = sapply(rq1x, function(var) {
      sum(!is.na(df[[var]]))
    }),
    `Study Wave` = sapply(rq1x, get_study_wave)
  )

v_rq1x$Description = str_remove(v_rq1x$Description,"\\(1st.*")
v_rq1x$Description = clean_description(v_rq1x$Description)

# Manual edit for cens01pop98density study wave
v_rq1x$`Study.Wave`[v_rq1x$`Teds.Code` == "cens01pop98density"] = "1st Contact"
v_rq1x$`Study.Wave`[v_rq1x$`Teds.Code` == "pollution1998pca"] = "1st Contact"

# Clean uniform descriptions for participation indicator variables
rq1y_twin_descriptions = c(
  "Parent-report twin booklet data present",
  "Web test data present",
  "Self-report questionnaire data present",
  "Self-report behaviour booklet data present",
  "Self-report questionnaire data present",
  "TEDS21 phase 1 self-report questionnaire data present",
  "TEDS26 Mental Health Questionnaire data present",
  "CATSLife web test data present"
)

v_rq1y =
  data.frame(
    `Short Label` = rq1y_twin_labels_clean,
    Description = rq1y_twin_descriptions,
    `Teds Code` = ifelse(rq1y_twin1 %in% colnames(df0), rq1y_twin1, paste0(rq1y_twin1, "*")),
    `Range or Level` = sapply(rq1y_twin1, function(var) {
      if (class(df[[var]]) == "numeric") {
        paste0(round(min(df[[var]], na.rm = TRUE), 2), " — ", round(max(df[[var]], na.rm = TRUE), 2))
      } else if (is.factor(df[[var]])) {
        factor_levels = levels(df[[var]])
        paste(c(paste0(factor_levels[1], "*"), factor_levels[-1]), collapse = ", ")
      } else {
        paste(unique(df[[var]]), collapse = ", ")
      }
    }),
    N = sapply(rq1y_twin1, function(var) {
      sum(!is.na(df[[var]]))
    }),
    `Study Wave` = sapply(rq1y_twin, get_study_wave)
  )

v_rq2 = data.frame(
  `Short Label` = rq2y_labels_short,
  Description   = rq2y_labels,
  `Teds Code` = ifelse(rq2y %in% colnames(df0), rq2y,paste0(rq2y,"*")),
  `Range or Level` = sapply(rq2y, function(var) {
    if (class(df[[var]]) == "numeric") {
      paste0(round(min(df[[var]], na.rm = TRUE), 2), " — ", round(max(df[[var]], na.rm = TRUE), 2))
    } else if (is.factor(df[[var]])) {
      factor_levels = levels(df[[var]])
      paste(c(paste0(factor_levels[1],"*"), factor_levels[-1]), collapse = ", ")
    } else {
      paste(unique(df[[var]]), collapse = ", ")
    }
  }),
  N = sapply(rq2y, function(var) {
    sum(!is.na(df[[var]]))
  }),
  `Study Wave` = sapply(rq2y, get_study_wave)
)

v_rq2$Description = str_remove(v_rq2$Description,"\\(2.*")

v_rq3 = data.frame(
  `Short Label` = rq6y_labels,
  Description   = clean_description(var_to_label(rq6y)),
  `Teds Code` = ifelse(rq6y %in% colnames(df0), rq6y, paste0(rq6y, "*")),
  `Range or Level` = sapply(rq6y, function(var) {
    if (class(df[[var]]) == "numeric") {
      paste0(round(min(df[[var]], na.rm = TRUE), 2), " — ", round(max(df[[var]], na.rm = TRUE), 2))
    } else if (is.factor(df[[var]])) {
      factor_levels = levels(df[[var]])
      paste(c(paste0(factor_levels[1], "*"), factor_levels[-1]), collapse = ", ")
    } else {
      paste(unique(df[[var]]), collapse = ", ")
    }
  }),
  N = sapply(rq6y, function(var) {
    sum(!is.na(df[[var]]))
  }),
  `Study Wave` = sapply(rq6y, get_study_wave)
)

v_rq5 = data.frame(
  `Short Label` = rq5y_labels_short,
  Description   = clean_description(var_to_label(rq5y)),
  `Teds Code` = ifelse(rq5y %in% colnames(df0), rq5y, paste0(rq5y, "*")),
  `Range or Level` = sapply(rq5y, function(var) {
    if (class(df[[var]]) == "numeric") {
      paste0(round(min(df[[var]], na.rm = TRUE), 2), " — ", round(max(df[[var]], na.rm = TRUE), 2))
    } else if (is.factor(df[[var]])) {
      factor_levels = levels(df[[var]])
      paste(c(paste0(factor_levels[1],"*"), factor_levels[-1]), collapse = ", ")
    } else {
      paste(unique(df[[var]]), collapse = ", ")
    }
  }),
  N = sapply(rq5y, function(var) {
    sum(!is.na(df[[var]]))
  }),
  `Study Wave` = sapply(rq5y, get_study_wave)
)

# Combine all v_rq dataframes with row indices for grouping
v_rq1x_indexed = cbind(row_group = "RQ1 Variables", row_id = 1:nrow(v_rq1x), v_rq1x)
v_rq1y_indexed = cbind(row_group = "RQ1 Outcome Variables", row_id = (nrow(v_rq1x)+1):(nrow(v_rq1x)+nrow(v_rq1y)), v_rq1y)
v_rq2_indexed = cbind(row_group = "RQ2 Variables", row_id = (nrow(v_rq1x)+nrow(v_rq1y)+1):(nrow(v_rq1x)+nrow(v_rq1y)+nrow(v_rq2)), v_rq2)
v_rq3_indexed = cbind(row_group = "RQ3 Variables", row_id = (nrow(v_rq1x)+nrow(v_rq1y)+nrow(v_rq2)+1):(nrow(v_rq1x)+nrow(v_rq1y)+nrow(v_rq2)+nrow(v_rq3)), v_rq3)
v_rq5_indexed = cbind(row_group = "RQ5 Variables", row_id = (nrow(v_rq1x)+nrow(v_rq1y)+nrow(v_rq2)+nrow(v_rq3)+1):(nrow(v_rq1x)+nrow(v_rq1y)+nrow(v_rq2)+nrow(v_rq3)+nrow(v_rq5)), v_rq5)

v_rq_combined = rbind(v_rq1x_indexed, v_rq1y_indexed, v_rq2_indexed, v_rq3_indexed, v_rq5_indexed)

gt(v_rq_combined) %>%
  tab_row_group(
    label = "E) RQ4 Outcome Variables",
    rows = row_group == "RQ5 Variables"
  ) %>%
  tab_row_group(
    label = "D) RQ3 Outcome Variables",
    rows = row_group == "RQ3 Variables"
  ) %>%
  tab_row_group(
    label = "C) RQ2 Outcome Variables",
    rows = row_group == "RQ2 Variables"
  ) %>%
  tab_row_group(
    label = "B) RQ1-2 Participation Indicators (1 = yes, 0 = no)",
    rows = row_group == "RQ1 Outcome Variables"
  ) %>%
  tab_row_group(
    label = "A) Baseline variables used to predict later participation, and for IP weighting",
    rows = row_group == "RQ1 Variables"
  ) %>%
  cols_hide(c(row_group, row_id)) %>%
  cols_width(
    `Range.or.Level` ~ px(100)
  ) %>%
  cols_label_with(fn = ~ gsub("\\.", " ", .x)) %>%
  tab_style(
    style = cell_text(font = "Times New Roman", size = px(12)),
    locations = cells_body()
  ) %>%
  tab_style(
    style = cell_text(font = "Times New Roman", size = px(12)),
    locations = cells_column_labels()
  ) %>%
  tab_style(
    style = cell_text(font = "Times New Roman", size = px(12)),
    locations = cells_row_groups()
  ) %>%
  tab_footnote(
    footnote = "Note. MFQ = Mood and Feelings Questionnaire; SDQ = Strengths and Difficulties Questionnaire; GCSE = General Certificate of Secondary Education; KS3 = Key Stage 3; GAD-D = Generalized Anxiety Disorder scale; PARCA = Parent Report of Children's Abilities; MZ = Monozygotic; DZ = Dizygotic; HNC = Higher National Certificate; HND = Higher National Diploma; CSE = Certificate of Secondary Education."
  )

Short Label	Description	Teds Code	Range or Level	N	Study Wave
A) Baseline variables used to predict later participation, and for IP weighting
Twin Sex	Twin sex, 0female 1male	sex1	0 — 1	26040	NA
Mother age at birth	Age in years of natural mother at time of birth of twins	amumagetw	14 — 45	25563	1st Contact
Father age at birth	Age in years of natural father at time of birth of twins	adadagetw	16 — 66	22975	1st Contact
Single Parent	Single Parent	asingle	cohabiting biological mother and father / cohabiting biological parent with other*, single parent	25587	1st Contact
Zygosity	Twin pair zygosity (best available estimate), 1MZ 2DZ	zygos	1*, 2	26040	26 Year
Mother medical risk	Mother medical risk factor composite scale	amedtot	-2.49 — 4.98	25864	1st Contact
Father employment level	Father Standard Occupational Classification employment level	afasoc2*	1*, 2, 3, 4, 5, 6, 7, 8, 9, caring for children at home, no job	23111	1st Contact
Father education level	Male parent highest qualification level	afahqual	CSE grade 1 or O-level/GCSE grade A-C*, no qualifications, CSE grade 2-5 or O-level/GCSE grade D-G, A-level or S-level, HNC, HND, undergraduate degree, postgraduate qualification	23049	1st Contact
Mother employment level	Mother Standard Occupational Classification employment level	amosoc2*	caring for children at home*, 1, 2, 3, 4, 5, 6, 7, 8, 9, no job	25650	1st Contact
Mother education level	Female parent highest qualification level	amohqual	no qualifications*, CSE grade 2-5 or O-level/GCSE grade D-G, CSE grade 1 or O-level/GCSE grade A-C, A-level or S-level, HNC, HND, undergraduate degree, postgraduate qualification	25629	1st Contact
Twin medical risk	Twin medical risk factor composite scale	atwmed1	-1.88 — 3.33	25832	1st Contact
Ethnic origin	Ethnic origin of twins, original codes	aethnicc	White*, Asian, Black, Mixed race, Other	25946	1st Contact
Language at home	Main language spoken at home	alang	other*, English, English + other	25678	1st Contact
Older siblings	Number of older siblings	anoldsib	0*, 1, 2, 3, 4, 5 or more	26040	1st Contact
Younger siblings	Number of younger siblings	anyngsib	0*, 1, 2 or more	26040	1st Contact
Twins club member	Member of a Twins Club	atwclub	0*, 1	25148	1st Contact
Childcare by others	Twins looked after by anyone else	alookels	0*, 1	24791	1st Contact
Smoking in pregnancy	Smoked cigarettes while pregnant	asmoke	0*, 1	25932	1st Contact
Alcohol in pregnancy	Drank alcohol while pregnant	adrink	0*, 1	25783	1st Contact
Severe stress in pregnancy	Severe stress during pregnancy	astress	0*, 1	25889	1st Contact
Pollution index	Principal Component of 1998 pollution variables	pollution1998pca	-2.47 — 4.57	23030	1st Contact
B) RQ1-2 Participation Indicators (1 = yes, 0 = no)
Y4 (parent-report twin booklet)	Parent-report twin booklet data present	dtwdata1	0 — 1	26040	4 Year
Y12 (web tests)	Web test data present	lcwdata1	0 — 1	26040	12 Year
Y12 (questionnaire)	Self-report questionnaire data present	lcqdata1	0 — 1	26040	12 Year
Y16 (behaviour booklet)	Self-report behaviour booklet data present	pcbhdata1	0 — 1	26040	16 Year
Y18 (questionnaire)	Self-report questionnaire data present	rcqdata1	0 — 1	26040	18 Year
Y21 (TEDS21 phase-1 questionnaire)	TEDS21 phase 1 self-report questionnaire data present	u1cdata1	0 — 1	26040	21 Year
Y26 (TEDS26 questionnaire)	TEDS26 Mental Health Questionnaire data present	zmhdata1	0 — 1	26040	26 Year
Y26 (CATSLife web tests)	CATSLife web test data present	zcdata1	0 — 1	26040	26 Year
C) RQ2 Outcome Variables
Maternal Education	Maternal Education (formatted as numeric variable)	amohqualn1	1 — 8	25629	1st Contact
Vocabulary	Vocabulary total score	bvocab1	0 — 100	11830	2 Year
Grammar	Grammar composite score	bgramma1	0 — 2	11783	2 Year
Parent-admin cognition	Parent-administered Parca mean score	badparn1	-3.09 — 3.45	11762	2 Year
Parent-report cognition	Parent-reported Parca total score	breparc1	0 — 26	11872	2 Year
Conduct problems	Conduct SDQ-comparable Behar subscale	bsdqccont1	0 — 8	11815	2 Year
Emotional problems	Emotion SDQ-comparable Behar subscale	bsdqcemot1	0 — 4	11823	2 Year
Hyperactivity	Hyperactivity SDQ-comparable Behar subscale	bsdqchypt1	0 — 6	11836	2 Year
Peer problems	Peer SDQ-comparable Behar subscale	bsdqcpert1	0 — 6	11565	2 Year
Prosocial behavior	Prosocial SDQ-comparable Behar subscale	bsdqcprot1	0 — 10	11729	2 Year
D) RQ3 Outcome Variables
Y12: Depression (MFQ)	MFQ scale from 11 MFQ items (child self-report) at 12	lcmfqt1	0 — 22	11432	12 Year
Y12: Externalising	SDQ Externalising scale at 12	lsdqext1	0 — 20	11389	12 Year
Y12: Cognitive ability	G composite scale from child web tests at 12, standardised	lcg1	-3.67 — 3.04	8458	12 Year
Y16: GCSE core subjects grade	Core subjects (English, maths, science): mean grade in GCSE results (twin exams at 16)	pcexgcsecoregrdm1	4 — 11	12982	16 Year
E) RQ4 Outcome Variables
Y12: Cognitive ability	G composite scale from child web tests at 12, standardised	lcg1	-3.67 — 3.04	8458	12 Year
Y14: Cognitive ability	G composite scale from child web tests at 14, standardised	ncg1	-4.12 — 3.17	5341	14 Year
Y16: Cognitive ability	G composite scale from child web tests at 16, standardised	pcg1	-2.86 — 4.06	4767	16 Year
Y21: G-game total score	G-game (general cognitive ability) overall total score	ucgt1	3 — 40	4549	21 Year
Y14: KS3 academic achievement	End of KS3 all-subject Academic achievement mean level (from parent questionnaire)	npks3tall1	1 — 9	5436	14 Year
Y16: GCSE core subjects grade	Core subjects (English, maths, science): mean grade in GCSE results (twin exams at 16)	pcexgcsecoregrdm1	4 — 11	12982	16 Year
Y21: Highest qualification	Twin probable highest level of qualification including current study (TEDS21 phase 1 twin questionnaire), 1-11	u1chqualp1	1 — 11	8912	21 Year
Y26: Highest qualification	Demographics item: highest qualification ordinal level (TEDS26 twin Mental Health Questionnaire)	zmhhqual1	1 — 11	8243	26 Year
Y12: Depression (MFQ)	MFQ scale from 11 MFQ items (child self-report) at 12	lcmfqt1	0 — 22	11432	12 Year
Y16: Depression (MFQ)	MFQ total scale (child behaviour questionnaire at 16)	pcbhmfqt1	0 — 26	9906	16 Year
Y21: Depression (MFQ)	MFQ overall total score (TEDS21 phase 1 twin questionnaire)	u1cmfqt1	0 — 16	9204	21 Year
Y26: Depression (MFQ)	MFQ overall total score (TEDS26 twin Mental Health Questionnaire)	zmhmfqt1	0 — 26	8306	26 Year
Y21: Anxiety (GAD-D)	General Anxiety overall total score (TEDS21 phase 2 twin questionnaire)	u2cganxt1	0 — 40	8236	21 Year
Y26: Anxiety (GAD-D)	GAD-D (General Anxiety) overall total score (TEDS26 twin Mental Health Questionnaire)	zmhganxt1	0 — 40	8022	26 Year
Y12: Externalising	SDQ Externalising scale at 12	lsdqext1	0 — 20	11389	12 Year
Y16: Externalising	SDQ Externalising scale at 16	psdqext1	0 — 20	9889	16 Year
Y21: Externalising	SDQ Externalising scale at 21	usdqext1	0 — 19	9210	21 Year
Y26: Externalising	SDQ Externalising scale at 26	zsdqext1	0 — 19	7718	26 Year
Note. MFQ = Mood and Feelings Questionnaire; SDQ = Strengths and Difficulties Questionnaire; GCSE = General Certificate of Secondary Education; KS3 = Key Stage 3; GAD-D = Generalized Anxiety Disorder scale; PARCA = Parent Report of Children's Abilities; MZ = Monozygotic; DZ = Dizygotic; HNC = Higher National Certificate; HND = Higher National Diploma; CSE = Certificate of Secondary Education.

Participation-Related plots

Missing Data at each timepoint

Family-Level Participation Metrics

df %>%
  select(c(acontact,rq1y)) %>%
  sapply(., function(x) length(which(x==1))/length(which(x>=0))) %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  `colnames<-`(c("Time Point","Participation Rate")) %>%
  mutate(`Time Point` = factor(`Time Point`, levels = `Time Point`, labels = var_to_label(`Time Point`))) %>%
  ggplot(aes(x = `Participation Rate`, y = `Time Point`)) + 
  geom_col() +
    geom_text(
    aes(label = paste0(round(`Participation Rate` * 100, 1), "%")),
    hjust = 1.1,  # Right align inside the bar (>1 moves it inside)
    color = "white",
    size = 3.5,
    fontface = "bold"
  ) +
  # geom_vline(aes(xintercept=.2)) +
  theme_bw()

Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
ℹ Please use `all_of()` or `any_of()` instead.
  # Was:
  data %>% select(rq1y)

  # Now:
  data %>% select(all_of(rq1y))

See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.

save_plot("11_participation_rate", width = 9, height = 5)


df %>%
  select(cohort, acontact, all_of(rq1y)) %>%
  group_by(cohort) %>%
  summarise(across(everything(), 
                   ~ length(which(. == 1)) / length(which(. >= 0)),
                   .names = "{.col}"),
            .groups = "keep") %>%
  pivot_longer(cols = -cohort, names_to = "Time Point", values_to = "Participation Rate") %>%
  mutate(
    `Time Point` = factor(`Time Point`, 
                              levels = c("acontact", rq1y), 
                              labels = var_to_label(c("acontact", rq1y))),
         Time_numeric = as.numeric(`Time Point`),
    cohort = gsub(" twins", "",cohort),
    cohort = gsub(" to ", "-",cohort)
    ) %>%
  ggplot(aes(x = Time_numeric, y = `Participation Rate`, color = cohort, group = cohort)) + 
  geom_point(size = 3) +
  geom_line(linewidth = 1) +
  geom_text(
    aes(label = paste0(round(`Participation Rate` * 100, 1), "%")),
    vjust = -0.5,
    hjust = 0.5,
    size = 3,
    show.legend = FALSE
  ) +
  scale_x_continuous(
    breaks = 1:length(c("acontact", rq1y)),
    labels = var_to_label(c("acontact", rq1y))
  ) +
  scale_y_continuous(labels = scales::percent_format()) +
  theme_bw() +
  labs(x = "Time Point", 
       y = "Participation Rate", 
       color = "Cohort",
       title = "Participation Rates Over Time by Cohort") +
  theme(legend.position = c(.8,.8),
        axis.text.x = element_text(angle = 45, hjust = 1)) + 
  theme(plot.margin = margin(t = 1, r = 1, b = 1, l = 10, unit = "mm"))

Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.

save_plot("11_participation_rate_by_cohort", width = 11, height = 9)


df %>%
  select(cohort, acontact, all_of(rq1y)) %>%
  filter(btwoyear==1) %>%
  group_by(cohort) %>%
  summarise(across(everything(), 
                   ~ length(which(. == 1)) / length(which(. >= 0)),
                   .names = "{.col}"),
            .groups = "keep") %>%
  pivot_longer(cols = -cohort, names_to = "Time Point", values_to = "Participation Rate") %>%
  mutate(
    `Time Point` = factor(`Time Point`, 
                              levels = c("acontact", rq1y), 
                              labels = var_to_label(c("acontact", rq1y))),
         Time_numeric = as.numeric(`Time Point`),
    cohort = gsub(" twins", "",cohort),
    cohort = gsub(" to ", "-",cohort)
    ) %>%
  ggplot(aes(x = Time_numeric, y = `Participation Rate`, color = cohort, group = cohort)) + 
  geom_point(size = 3) +
  geom_line(linewidth = 1) +
  geom_text(
    aes(label = paste0(round(`Participation Rate` * 100, 1), "%")),
    vjust = -0.5,
    hjust = 0.5,
    size = 3,
    show.legend = FALSE
  ) +
  scale_x_continuous(
    breaks = 1:length(c("acontact", rq1y)),
    labels = var_to_label(c("acontact", rq1y))
  ) +
  scale_y_continuous(labels = scales::percent_format()) +
  theme_bw() +
  labs(x = "Time Point", 
       y = "Participation Rate", 
       color = "Cohort",
       title = "Participation Rates Over Time by Cohort",
       subtitle = "Just selecting participants that participated at YEAR 2"
       ) +
  theme(legend.position = c(.8,.8),
        axis.text.x = element_text(angle = 45, hjust = 1)) + 
  theme(plot.margin = margin(t = 1, r = 1, b = 1, l = 10, unit = "mm"))

save_plot("11_participation_rate_by_cohort_selectedy2", width = 11, height = 9)

Twin-Level Participation Metrics

Participation Rates

df %>%
  select(cohort, all_of(rq1y_twin1)) %>%
  # select(cohort, !starts_with(vars_to_na)) %>%
  group_by(cohort) %>%
  summarise(across(everything(), 
                   ~ length(which(. == 1)) / length(which(. >= 0)),
                   .names = "{.col}"),
            .groups = "keep") %>%
  pivot_longer(cols = -cohort, names_to = "Time Point", values_to = "Participation Rate") %>%
  mutate(
    `Time Point` = factor(`Time Point`, 
                              levels = paste0(rq1y_twin,"1"), 
                              labels = rq1y_twin_labels),
    # Time_numeric = as.numeric(`Time Point`),
    cohort = gsub(" twins", "",cohort),
    cohort = gsub(" to ", "-",cohort)
    ) %>%
  ggplot(aes(x = `Time Point`, y = `Participation Rate`, color = cohort, group = cohort)) + 
  geom_point(size = 3) +
  geom_line(linewidth = 1) +
  geom_text(
    aes(label = paste0(round(`Participation Rate` * 100, 1), "%")),
    vjust = -0.5,
    hjust = 0.5,
    size = 3,
    show.legend = FALSE
  ) +
  # scale_x_continuous(
  #   breaks = 1:length(c("acontact", rq1y)),
  #   labels = var_to_label(c("acontact", rq1y))
  # ) +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0,1)) +
  theme_bw() +
  labs(x = "Time Point", 
       y = "Participation Rate", 
       color = "Cohort",
       title = "Participation Rates Over Time by Cohort") +
  theme(legend.position = c(.8,.8),
        axis.text.x = element_text(angle = 45, hjust = 1)) + 
  theme(plot.margin = margin(t = 1, r = 1, b = 1, l = 10, unit = "mm")) 

Twin Concordance

NOTE THAT WE REMOVE COHORTS 2 & 3

ALSO FOR CORRELATIONS - I SHOULD FILTER BY TWIN == 1!

vars_to_na = c("btwdata","ctwdata", "gcdata","icdata","jcdata","pcwebdata")

df_cohort34_removed = df %>%
  mutate(
      across(starts_with(vars_to_na), ~ ifelse(cohort %in% c("Cohort 3: twins born Sep-95 to Aug-96","Cohort 4: twins born Sep-96 to Dec-96"), NA, .))
    )

plot_data = map_dfr(rq1y_twin, function(var){
  df_cohort34_removed %>% 
    group_by(sexzyg) %>%
    summarise(
      mean = mean(.data[[paste0(var,"1")]], na.rm = TRUE),
      cor  = cor(x=.data[[paste0(var,"1")]], y=.data[[paste0(var,"2")]], use = "pairwise.complete.obs")
    # cor  = {
    #   group_data <- cur_data()
    #   as.numeric(psych::tetrachoric(x=group_data[paste0(var,1:2)])$rho)[2]
    # }
    ) %>%
    ungroup() %>%
    mutate(variable = var)
}, .id = "var_index") %>%
  pivot_longer(cols = c(mean, cor), 
               names_to = "metric", 
               values_to = "value") %>%
  mutate(
    metric_label = case_when(
      metric == "mean" ~ "Participation Rate",
      metric == "cor" ~ "Twin Tetrachoric Correlation"
    ),
    time_point = rq1y_twin_labels_clean[match(.$variable, rq1y_twin)],
    time_point = factor(time_point, levels = rq1y_twin_labels_clean)
  )

plot_data %>%
  ggplot( aes(x = time_point, y = value, color = sexzyg, group = sexzyg)) +
  geom_line(size = 1, alpha = 0.8) +
  geom_point(size = 2.5, alpha = 0.9) +
  ggrepel::geom_text_repel(aes(label = case_when(
    metric == "mean" ~ paste0(round(value * 100, 1), "%"),
    metric == "cor" ~ gbtoolbox::apa_num(value)
  )),
  position = position_dodge(width = 0.5),
  vjust = 2,
  size = 2,
  show.legend = FALSE) +
  facet_wrap(~metric_label, scales = "free_y", ncol = 1) +
  scale_color_manual(
    values = c(
      "MZ male" = "#d7191c",
      "DZ male" = "#fdae61",
      "MZ female" = "#2c7bb6",
      "DZ female" = "#abd9e9",
      "DZ opposite sexes" = "grey"
    )
  ) +
  ggh4x::scale_y_facet(metric_label == "Participation Rate",
                       labels = scales::percent_format(),
                       breaks = seq(0, 1, 0.1)) +
  labs(
    # title = "Twin Study: Participation Rates and Correlations Across Time Points",
    # subtitle = "Trends by Zygosity and Sex",
    x = "Participation Time Point",
    y = NULL,
    color = "Zygosity",
  ) +
  theme_minimal() +
  theme(
    # plot.title = element_text(size = 14, face = "bold"),
    # plot.subtitle = element_text(size = 12),
    # strip.text = element_text(size = 11, face = "bold"),
    legend.spacing.y = unit(0, "pt"),
    legend.background = element_rect(fill = alpha("white",.5), color = NA),
    legend.key.height = unit(0.6, "pt"),
    legend.text = element_text(margin = margin(t = 0, b = 0)),
    legend.position = c(0.02, 0.02),
    legend.justification = c(0, 0),
    # legend.title = element_text(face = "bold"),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
    plot.margin = margin(t = 5.5, r = 5.5, b = 5.5, l = 47, unit = "pt")
  ) +
  guides(color = guide_legend(override.aes = list(size = 3)))

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

save_plot("11_participation_rates_zygosity", width = 7, height = 6)

Create ACE estimates plot

estimates = readRDS(file.path("results","1_ace_estimates.Rds"))
ci       = readRDS(file.path("results","1_ace_confints.Rds"))

estimates %>%
    mutate(
      par_clean = case_when(
        par == "a1" ~ "A",
        par == "c1" ~ "C",
        par == "e1" ~ "E"
      ),
      time_point = rq1y_twin_labels_clean[match(name, rq1y_twin)],
      time_point = factor(time_point, levels = rq1y_twin_labels_clean)
    ) %>%
    ggplot(aes(x = time_point, y = value, color = par_clean, group = interaction(sex, par_clean))) +
    geom_line(size = .5, alpha = 0.8, aes(linetype = sex)) +
    geom_point(size = 2.5, alpha = 0.9, aes(shape = sex)) +
    # ggrepel::geom_text_repel(aes(label = sprintf("%.2f", value)),
    #   size = 1.5,
    #   show.legend = FALSE,
    #   min.segment.length = 0,
    #   segment.size = 0.2,
    #   max.overlaps = Inf,
    #   direction = "y",
    #   force = 2) +
    scale_color_manual(
      values = c(
        "A" = "#d7191c",
        "C" = "#fdae61", 
        "E" = "#2c7bb6"
      )
    ) +
    scale_linetype_manual(
      values = c("male" = "solid", "female" = "dotted")
    ) +
    scale_shape_manual(
      values = c("male" = 16, "female" = 17)
    ) +
    scale_y_continuous(labels = scales::number_format(accuracy = 0.01)) +
    labs(
      x = "Participation Time Point",
      y = NULL,
      color = "ACE Component",
      linetype = "Sex",
      shape = "Sex"
    ) +
    theme_minimal() +
    theme(
      legend.background = element_blank(),
      legend.position = "top",
      legend.justification = "center",
      legend.direction = "horizontal",
      legend.box = "horizontal",
      legend.key.height = unit(0.5, "lines"),
      legend.key.width = unit(0.8, "lines"),
      legend.text = element_text(size = 8),
      legend.title = element_text(size = 9),
      legend.margin = margin(2, 2, 2, 2),
      panel.grid.minor = element_blank(),
      axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
      plot.margin = margin(t = 5.5, r = 5.5, b = 5.5, l = 47, unit = "pt"),
      plot.title = element_text(size = 8, hjust = 0.5),
      strip.text = element_text(size = 10, face = "bold")
    ) +
    guides(
      color = guide_legend(override.aes = list(size = 2), order = 1),
      linetype = guide_legend(override.aes = list(size = 0.8), order = 2),
      shape = guide_legend(override.aes = list(size = 1.5), order = 2)
    )

Warning: Duplicated `override.aes` is ignored.

  save_plot("11_ACE_time", width = 7, height = 3)

Warning: Duplicated `override.aes` is ignored.
Duplicated `override.aes` is ignored.

ANOVA Sex & Zygosity predicting participation

map_dfr(seq_along(rq1y_twin1), function(i){
  y_var     = rq1y_twin1[i]
  model     = glm(paste(y_var, "~ zygos + sex1"), data = df, family = binomial)
  robust_se = clubSandwich::coef_test(model, cluster = df$randomfamid, vcov = "CR0", test = "z")
  
  data.frame(
    `Time Point`    = rq1y_twin_labels_clean[i],
    `Zygosity Coef` = as.numeric(robust_se[2, "beta"]),
    `Zygosity p`    = as.numeric(robust_se[2, "p_z"]),
    `Sex Coef`      = as.numeric(robust_se[3, "beta"]),
    `Sex p`         = as.numeric(robust_se[3, "p_z"])
  )
}) %>%
  knitr::kable(digits = 3)

Registered S3 method overwritten by 'clubSandwich':
  method    from    
  bread.mlm sandwich

Time.Point	Zygosity.Coef	Zygosity.p	Sex.Coef	Sex.p
Y4 (parent-report twin booklet)	-0.064	0.093	-0.098	0.001
Y12 (web tests)	-0.240	0.000	-0.241	0.000
Y12 (questionnaire)	-0.225	0.000	-0.176	0.000
Y16 (behaviour booklet)	-0.181	0.000	-0.338	0.000
Y18 (questionnaire)	-0.169	0.000	-0.224	0.000
Y21 (TEDS21 phase-1 questionnaire)	-0.201	0.000	-0.752	0.000
Y26 (TEDS26 questionnaire)	-0.193	0.000	-0.875	0.000
Y26 (CATSLife web tests)	-0.258	0.000	-0.790	0.000

How many different time points do we have for each participant (excluding first contact)

df %>% 
  select(c(rq1y)) %>% 
  apply(., 1, function(x) sum(x)) %>% 
  tibble(`Number Time Points Participated` = .) %>% 
  count(`Number Time Points Participated`) %>% 
  ggplot(aes(x = `Number Time Points Participated`, y = n)) + 
  geom_col() + 
  geom_text(aes(label = n), vjust = 1.2, color = "white", size = 3.5) +
  theme_bw() + 
  labs(
    title = "How many times has each pps been tested (exl. first contact)?",
    subtitle = "Twin-level participation numbers (N = 27890)"
    )

save_plot("11_number_of_times_each_twin_participated", width = 9, height = 5)

How predictive is cohort and year of birth in participation at each timepoint?

I’ve checked there are no missing values for cohort or year of birth (aonsby)

View full-sizeplot

models = list()

for (i in seq_along(rq1y)){
  # cat(i, "/", length(rq1y),"\n")
  formula <- as.formula(paste(rq1y[i], "~ cohort + aonsby"))
  models[[i]] = glm(formula, data = df, family = binomial)
}

var_explained = sapply(models, function(x) as.numeric(performance::r2(x)))

cbind.data.frame(rq1y, rq1y_labels, var_explained) %>%
  `colnames<-`(c("Var","Label","R2 by Cohort and Birth Year")) %>%
  knitr::kable(digits = 2)

Var	Label	R2 by Cohort and Birth Year
btwoyear	2 Year data present, 1Y 0N	0.37
cthreeyr	3 Year data present, 1Y 0N	0.34
dfouryr	4 Year data present, 1Y 0N	0.01
gsevenyr	7 year study data present, 1Y 0N	0.00
heightyr	8yr data present, 1Y 0N	0.01
jtenyear	10 year study data present, 1Y 0N	0.31
ltwelvyr	12yr data present, 1Y 0N	0.01
n14year	14yr data present, 1Y 0N	0.02
p16year	16 Year study data present, 1Y 0N	0.01
rcqdata	18 Year questionnaire data present (either twin), 1Y 0N	0.05
uteds21data	TEDS21 study data present (parent or either twin, either phase), 1Y 0N	0.00
zmhdata	TEDS26 twin MHQ data present for at least one twin in the pair, 1Y 0N	0.00

cbind.data.frame(rq1y, rq1y_labels, var_explained) %>%
  `colnames<-`(c("Var","Label","Variance Explained")) %>%
  mutate(
    Label = factor(Label, levels = rq1y_labels)
  ) %>%
  ggplot(aes(x=Label,y=`Variance Explained`)) + 
  geom_col() + 
  geom_text(aes(label = gbtoolbox::apa_num(`Variance Explained`)),
            vjust = 1.1, color = "white", size = 3) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust=1, vjust=1)) +
  scale_x_discrete(expand = expansion(mult = c(0.1, 0)))

save_plot("11_r2_participation_cohort_birthyear", width = 9, height = 5)

At the Y18 time point, we see that those in cohort 1 and those born in 1994 are less likely to participate.

rq1y[10]

[1] "rcqdata"

var_to_label(rq1y[10])

[1] "18 Year questionnaire data present (either twin), 1Y 0N"

summary(models[[10]])

Call:
glm(formula = formula, family = binomial, data = df)

Coefficients:
                                            Estimate Std. Error z value
(Intercept)                                 -0.63412    0.02676 -23.696
cohortCohort 2: twins born Sep-94 to Aug-95  0.93646    0.04510  20.766
cohortCohort 3: twins born Sep-95 to Aug-96  0.87605    0.06398  13.692
cohortCohort 4: twins born Sep-96 to Dec-96  0.72585    0.07991   9.083
aonsby1995                                   0.18654    0.04495   4.150
aonsby1996                                   0.14914    0.06514   2.290
                                            Pr(>|z|)    
(Intercept)                                  < 2e-16 ***
cohortCohort 2: twins born Sep-94 to Aug-95  < 2e-16 ***
cohortCohort 3: twins born Sep-95 to Aug-96  < 2e-16 ***
cohortCohort 4: twins born Sep-96 to Dec-96  < 2e-16 ***
aonsby1995                                  3.32e-05 ***
aonsby1996                                     0.022 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 35954  on 26039  degrees of freedom
Residual deviance: 34729  on 26034  degrees of freedom
AIC: 34741

Number of Fisher Scoring iterations: 4

Cohort alone is not helpful for determining who was eligible at t2 and t3

As we can see below, nobody in cohort 4 took part in Y2, Y3, Y10 data collection. However, a minority of participants took part in Y2-Y3 from cohort 3, but participation rates are quite low suggesting that some participants were elibile and some not.

 table(df0$cohort, df0$btwoyear)

   
       0    1
  1 1866 4664
  2 3900 5876
  3 6558 1978
  4 3048    0

 table(df0$cohort, df0$cthreeyr)

   
       0    1
  1 2348 4182
  2 3838 5938
  3 6622 1914
  4 3048    0

 table(df0$cohort, df0$jtenyear)

   
       0    1
  1 3794 2736
  2 4390 5386
  3 8536    0
  4 3048    0

 table(df0$AONSBY, df0$btwoyear)

      
          0    1
  1994 2788 7036
  1995 4320 5482
  1996 8264    0

 table(df0$AONSBY, df0$cthreeyr)

      
          0    1
  1994 3678 6146
  1995 3914 5888
  1996 8264    0

 table(df0$AONSBY, df0$jtenyear)

      
          0    1
  1994 5288 4536
  1995 6216 3586
  1996 8264    0

Participation Patterns

df %>%
  select(c(rq1y)) %>% 
  rename_with(~var_to_label(., df)) %>%
  mutate_all(~na_if(., 0)) %>%
  ggmice::plot_pattern(npat = 70, rotate = TRUE) + 
  scale_x_discrete(limits = (var_to_label(rq1y))) + 
  coord_flip()

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Coordinate system already present. Adding new coordinate system, which will
replace the existing one.

save_plot("11_participation_pattern", width = 19, height = 6)


df %>%
  select(c(rq1y)) %>%
  mutate_all(~ .==1) %>%
  rowwise() %>%
  mutate(
    `Year 2 - 4 present`   = any(btwoyear, cthreeyr, dfouryr),
    `Year 7 - 12 present`  = any(gsevenyr, heightyr, jtenyear),
    `Year 14 - 16 present` = any(n14year, p16year),
    `Year 18 - 26 present` = any(rcqdata, uteds21data, zmhdata)
  ) %>%
  select(ends_with("present")) %>%
  mutate_all(~ na_if(., FALSE)) %>%
  ggmice::plot_pattern(rotate = TRUE) + 
  scale_x_discrete(limits = c("Year 2 - 4 present", "Year 7 - 12 present", "Year 14 - 16 present", "Year 18 - 26 present")) + 
  coord_flip()

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Coordinate system already present. Adding new coordinate system, which will
replace the existing one.

save_plot("11_participation_pattern_grouped", width = 19, height = 6)

df %>%
  select(c(rq1y)) %>%
  rename_with(~var_to_label(., df)) %>%
  rename_with(~clean_rq1y_label(.)) %>%
  mutate_all(~ na_if(., 0)) %>%
  gbtoolbox::plot_pairwise_missing(textadjust = 2.5) +
  labs(subtitle = "Number of twins tested at each time point (diagonal) and \npair of time points (below diagonal)")

Warning in gbtoolbox::plot_pairwise_missing(., textadjust = 2.5): This function is in early beta, and not yet ready for widespread use. 
  Proceed with caution

save_plot("11_pairwise_missingness_plot", width = 6, height = 6)

Descriptive stats on specific variables

Age of twins at each time point

View full-size plot

df_age = df %>%
  select(contains("age")) %>%
  select(-contains("genpro"), -contains("mumage"), - contains("dadage"), -ends_with("2")) 

df_age_long = df_age %>%
  `colnames<-`(var_to_label(colnames(df_age))) %>%
  pivot_longer(cols = everything())

df_age_long %>% 
  mutate(name = factor(name, levels = var_to_label(colnames(df_age)))) %>%
  ggplot(aes(x=value)) + 
  geom_histogram(bins = 100) +
  facet_wrap(~name, ncol =3) + 
  scale_x_continuous(breaks = seq(0,30,by=2)) +
  theme_bw() +
  theme(
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank()
    )

Warning: Removed 202363 rows containing non-finite outside the scale range
(`stat_bin()`).

save_plot("11_participant_ages_each_timepoint", width = 12, height = 12)

Warning: Removed 202363 rows containing non-finite outside the scale range
(`stat_bin()`).
Removed 202363 rows containing non-finite outside the scale range
(`stat_bin()`).

Education

all(table(df$amohqual)==table(df$amohqualn))

[1] TRUE

df %>%
  filter(twin ==1) %>%
  count(amohqual) %>%
  mutate(
    percent = n/sum(n)
  ) %>%
  gt() %>% 
  fmt_percent(columns = percent)

Female parent highest qualification level (1st Contact), see value labels	n	percent
no qualifications	1208	9.28%
CSE grade 2-5 or O-level/GCSE grade D-G	2000	15.36%
CSE grade 1 or O-level/GCSE grade A-C	4779	36.71%
A-level or S-level	1706	13.10%
HNC	359	2.76%
HND	486	3.73%
undergraduate degree	1495	11.48%
postgraduate qualification	782	6.01%
NA	205	1.57%

Relationship between maternal education and child outcomes

This section examines whether maternal education level predicts child academic and cognitive outcomes. We regress each outcome on maternal education using dummy coding (treating each education level as a separate category), which allows us to detect non-linear relationships.

# Helper function to create marginal effects plot
plot_education_effects = function(model, y_label) {
  marginaleffects::predictions(model, by = "amohqual") %>%
    data.frame() %>%
    rename(
      prediction = estimate,
      education_level = amohqual
    ) %>%
    ggplot(aes(y = prediction, x = education_level, group = 1)) +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
    geom_point(size = 2) +
    geom_line() +
    labs(
      x = "Maternal Education Level",
      y = y_label
    ) +
    theme_bw() +
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)
    )
}

KS3 Academic Achievement (Age 14)

model_ks3 = df %>%
  filter(twin == 1) %>%
  filter(!is.na(amohqual)) %>%
  lm(npks3tall1 ~ amohqual, data = .)

summary(model_ks3)

Call:
lm(formula = npks3tall1 ~ amohqual, data = .)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.8335 -0.4234  0.0387  0.4676  3.4857 

Coefficients:
                                                Estimate Std. Error t value
(Intercept)                                      5.23404    0.06752  77.516
amohqualCSE grade 2-5 or O-level/GCSE grade D-G  0.28023    0.08062   3.476
amohqualCSE grade 1 or O-level/GCSE grade A-C    0.53886    0.07145   7.541
amohqualA-level or S-level                       0.78125    0.07623  10.248
amohqualHNC                                      0.66686    0.10415   6.403
amohqualHND                                      0.66195    0.09346   7.082
amohqualundergraduate degree                     0.97547    0.07648  12.754
amohqualpostgraduate qualification               1.01003    0.08401  12.023
                                                Pr(>|t|)    
(Intercept)                                      < 2e-16 ***
amohqualCSE grade 2-5 or O-level/GCSE grade D-G 0.000517 ***
amohqualCSE grade 1 or O-level/GCSE grade A-C   6.32e-14 ***
amohqualA-level or S-level                       < 2e-16 ***
amohqualHNC                                     1.79e-10 ***
amohqualHND                                     1.80e-12 ***
amohqualundergraduate degree                     < 2e-16 ***
amohqualpostgraduate qualification               < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7397 on 2693 degrees of freedom
  (10114 observations deleted due to missingness)
Multiple R-squared:  0.1102,    Adjusted R-squared:  0.1078 
F-statistic: 47.63 on 7 and 2693 DF,  p-value: < 2.2e-16

plot_education_effects(model_ks3, "Predicted KS3 Score")

save_plot("11_meducation_ks3_plot", width = 6, height = 5)

GCSE Core Subjects Grade (Age 16)

model_gcse = df %>%
  filter(twin == 1) %>%
  filter(!is.na(amohqual)) %>%
  lm(pcexgcsecoregrdm1 ~ amohqual, data = .)

summary(model_gcse)

Call:
lm(formula = pcexgcsecoregrdm1 ~ amohqual, data = .)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.8174 -0.6731  0.0569  0.8087  3.0915 

Coefficients:
                                                Estimate Std. Error t value
(Intercept)                                      7.90847    0.06435 122.894
amohqualCSE grade 2-5 or O-level/GCSE grade D-G  0.28280    0.07571   3.735
amohqualCSE grade 1 or O-level/GCSE grade A-C    0.76463    0.06818  11.214
amohqualA-level or S-level                       1.20897    0.07350  16.449
amohqualHNC                                      1.13466    0.10065  11.274
amohqualHND                                      1.21039    0.09213  13.137
amohqualundergraduate degree                     1.73733    0.07321  23.729
amohqualpostgraduate qualification               1.79211    0.08073  22.198
                                                Pr(>|t|)    
(Intercept)                                      < 2e-16 ***
amohqualCSE grade 2-5 or O-level/GCSE grade D-G 0.000189 ***
amohqualCSE grade 1 or O-level/GCSE grade A-C    < 2e-16 ***
amohqualA-level or S-level                       < 2e-16 ***
amohqualHNC                                      < 2e-16 ***
amohqualHND                                      < 2e-16 ***
amohqualundergraduate degree                     < 2e-16 ***
amohqualpostgraduate qualification               < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.105 on 6431 degrees of freedom
  (6376 observations deleted due to missingness)
Multiple R-squared:  0.1821,    Adjusted R-squared:  0.1812 
F-statistic: 204.5 on 7 and 6431 DF,  p-value: < 2.2e-16

plot_education_effects(model_gcse, "Predicted GCSE Score")

save_plot("11_meducation_gcse_plot", width = 6, height = 5)

Cognitive Ability (Age 12)

model_cog = df %>%
  filter(twin == 1) %>%
  filter(!is.na(amohqual)) %>%
  lm(lcg1 ~ amohqual, data = .)

summary(model_cog)

Call:
lm(formula = lcg1 ~ amohqual, data = .)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8039 -0.6214  0.0563  0.6778  2.6524 

Coefficients:
                                                Estimate Std. Error t value
(Intercept)                                     -0.54041    0.06525  -8.282
amohqualCSE grade 2-5 or O-level/GCSE grade D-G  0.23722    0.07851   3.022
amohqualCSE grade 1 or O-level/GCSE grade A-C    0.43448    0.06979   6.225
amohqualA-level or S-level                       0.67126    0.07566   8.872
amohqualHNC                                      0.58609    0.10728   5.463
amohqualHND                                      0.65015    0.09575   6.790
amohqualundergraduate degree                     0.92674    0.07491  12.372
amohqualpostgraduate qualification               0.93773    0.08317  11.275
                                                Pr(>|t|)    
(Intercept)                                      < 2e-16 ***
amohqualCSE grade 2-5 or O-level/GCSE grade D-G  0.00253 ** 
amohqualCSE grade 1 or O-level/GCSE grade A-C   5.28e-10 ***
amohqualA-level or S-level                       < 2e-16 ***
amohqualHNC                                     4.95e-08 ***
amohqualHND                                     1.28e-11 ***
amohqualundergraduate degree                     < 2e-16 ***
amohqualpostgraduate qualification               < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9634 on 4195 degrees of freedom
  (8612 observations deleted due to missingness)
Multiple R-squared:  0.07031,   Adjusted R-squared:  0.06876 
F-statistic: 45.32 on 7 and 4195 DF,  p-value: < 2.2e-16

plot_education_effects(model_cog, "Predicted Cognitive Score (g)")

save_plot("11_meducation_cognitive_plot", width = 6, height = 5)